Optimization hardness in biological transients

Link to live notebook

Link to static webpage

Paper

Preamble

• After opening this notebook, run the cells below before anything else. These will import libraries and define functions.

• This notebook requires the following packages, which are pre-installed on Google Colab.

  • Python 3
  • NumPy
  • SciPy
  • Matplotlib

## Preamble / required packages
import numpy as np
from IPython.display import Image
np.random.seed(0)

## Import local plotting functions and in-notebook display functions
import matplotlib.pyplot as plt
from IPython.display import Image, display
%matplotlib inline
## default first plot color is black
plt.rcParams['axes.prop_cycle'] = plt.cycler(color=['k'])  # Set default plot color to black
plt.rcParams['axes.prop_cycle'] = plt.cycler(color=plt.cm.tab10.colors)  # Reset to default colors
plt.rcParams['legend.frameon'] = False

## Disable warnings
import warnings
warnings.filterwarnings('ignore')

# !apt-get install nvidia-cuda-toolkit
# !pip3 install numba
## These lines allow Google Colab to import and work with external files.
# %rm -rf illotka
# !git clone https://github.com/williamgilpin/illotka.git
# %cd illotka
# import importlib.util, sys
# module_path = "/content/illotka/base.py"
# spec = importlib.util.spec_from_file_location("base", module_path)
# base = importlib.util.module_from_spec(spec)
# sys.modules["base"] = base
# spec.loader.exec_module(base)

# import requests
# url = "https://raw.githubusercontent.com/williamgilpin/illotka/main/base.py"
# resp = requests.get(url)
# resp.raise_for_status()
# with open("base.py", "w") as f:
#     f.write(resp.text)

# url = "https://raw.githubusercontent.com/williamgilpin/illotka/main/utils.py"
# resp = requests.get(url)
# resp.raise_for_status()
# with open("utils.py", "w") as f:
#     f.write(resp.text)


from base import *
from utils import *

# Autoreload
%load_ext autoreload
%autoreload 2
print("Loading complete")

Loading complete

Revisiting stability vs. complexity

• Biological systems

• May 1972: We expect large ecosystems to be unstable. If we assume that the interaction matrix $A$ is drawn from a normal distribution, then the steady-state solution is given by

$$ \frac{dN_i}{dt} = \sum_{j=1}^N A_{ij} N_j $$

where $\mathbf{N} = (N_1, N_2, \cdots, N_N)$ is the vector of species abundances, and $A_{ij}$ is the interaction coefficient between species $i$ and $j$.

A = np.random.normal(size=(n, n))
eigvals = np.linalg.eigvals(A)

plt.figure(figsize=(5, 5))
plt.plot(eigvals.real, eigvals.imag, '.')

# np.max(np.real(np.linalg.eigvals(A)))

[<matplotlib.lines.Line2D at 0x4a8557ed0>]

The generalized Lotka-Volterra model

We consider random ecosystems given by the generalized Lotka-Volterra equation,

$$ \frac{dN_i}{dt} = N_i \left( r_i + \sum_{j=1}^N A_{ij} N_j \right) $$

where $N_i$ is the population of species $i$, $r_i$ is the intrinsic growth rate of species $i$, and $A_{ij}$ is the interaction coefficient between species $i$ and $j$.

To study the behavior of this model, we consider the case where $r_i \sim \mathcal{N}(0,1)$, and the matrix $A$ has the form

$$ A = Q - d\, I $$ where $I \in \mathbb{R}^{N \times N}$ is the identity matrix, $Q_{ij} \sim \mathcal{N}(0,1)$ and $d$ is a constant density-limitation.

# Set the simulation parameters
tmax = 20000 # Maximum integration duration
n = 200 # Number of species
d = 12.5 # Density-limitation

## Initialize the model
eq = GaussianLotkaVolterra(n)
eq.A = np.random.normal(size=(n, n)) - d * np.identity(n)
eq.r = np.random.normal(size=n)

## Initial conditions
ic = 0.4 *np.random.uniform(size=n)

## Numerically integrate the model
t, y = eq.integrate(tmax, ic)

plt.figure(figsize=(7, 4))
plt.semilogx(t, y, color="k", lw=1.5, alpha=0.3);
plt.xlabel("Time")
plt.ylabel("Species abundance")
plt.xlim(1e-2, np.max(t))
plt.ylim(0, None)

(0.0, 0.6341573011172086)

How many species survive?

nonzero_remaining = np.sum(y[-1] > 1e-12) # Precision floor
print(f"{nonzero_remaining} species out of {n} remain")

106 species out of 200 remain

We repeat the experiment many times

n_replicates = 20
tmax = 1000

number_survive = []
for i in prange(n_replicates):
    progress_bar(i, n_replicates) # Print progress bar for loop

    eq = GaussianLotkaVolterra(n)
    eq.A = np.random.normal(size=(n, n)) - d * np.identity(n)
    # eq.A = eq.A * np.identity(n) # Wipe out all off-diagonal elements
    eq.r = np.random.normal(size=n)
    t, y = eq.integrate(tmax, ic)

    n_survive = np.sum(y[-1] > 1e-12)
    number_survive.append(n_survive)

plt.hist(number_survive)
plt.xlabel("Number of surviving species");
plt.ylabel("Frequency");

[################### ] 

print(np.mean(number_survive))
print(np.std(number_survive))

103.66
6.259744403727679

print(0.5 * n)
print(np.sqrt(n * 0.5 * (1 - 0.5)))

100.0
7.0710678118654755

Properties of the Random Lotka-Volterra model

Proven in Serván et al. 2018.

1. On average, half of all species survive at steady-state. Across replicates, the number of surviving species follows a binomial distribution with mean $n/2$.
2. For a given matrix $A$, the steady-state solution is unique (global attractor).
3. For a given matrix $A$, there always exists a density-limitation $d$ that that leads to a single stable steady-state.

Main assumption: Interactions and growth rates are drawn from a symmetric distribution with finite variance.

How do extreme values affect the dynamics?

The generalized normal distribution is parametrized by shape parameter $\beta$, which controls the tail thickness. It is given by

$$ p(x) = \frac{\beta}{2\sigma\Gamma(1/\beta)} e^{-\left(\frac{|x-\mu|}{\sigma}\right)^\beta} $$

• when $\beta < 1$, the distribution has fatter tails than the normal distribution.
• when $\beta > 1$, the distribution has lighter tails than the normal distribution.

# Set the simulation parameters
tmax = 10000 # Maximum integration duration
n = 300 # Number of species
d = 18 # Density-limitation
seed = 0 # Seed for random number generator
# seed = None

## Initialize the model
np.random.seed(seed)
eq = GaussianLotkaVolterra(n, random_state=seed, tolerance=1e-10)

## Sample the interaction matrix from a normal distribution
# eq.A = np.random.normal(size=(n, n)) - d * np.identity(n)
# eq.r = np.random.normal(size=n)

## Sample the interaction matrix from a modified normal distribution
eq.A = 1.5 * normal_generalized(size=(n, n)) - d * np.identity(n)
eq.r = normal_generalized(size=n)

## Random initial conditions
ic = 0.5 * np.random.uniform(size=n)

## Numerically integrate the model
t, y = eq.integrate(tmax, ic)

## Plot the results
plt.figure(figsize=(7, 4))
plt.semilogx(t, y, color="k", lw=1.5, alpha=0.2);
plt.xlabel("Time")
plt.ylabel("Species abundance")
plt.xlim(1e-2, 1000)
plt.ylim(0, 0.1)


nonzero_remaining = np.sum(y[-1] > 1e-12) # Precision floor
print(f"{nonzero_remaining} species out of {n} remain")

156 species out of 300 remain

How long to reach steady-state?

Given a steady-state solution $\mathbf{N}^*$, the settling time is the time it takes for the solution to reach the steady-state within a given tolerance $\epsilon$.

$$ \text{argmin}_{t} \bigg( \dfrac{\| \mathbf{N}(t) - \mathbf{N}^* \|}{\| \mathbf{N}(t) \|} < \epsilon \bigg) $$

def compute_settling_time(traj, tvals=None, eps=1e-6):
    """
    Compute the settling time from a multivariate time series. 
    If the trajectory does not reach the steady-state within the integration time, 
    return a NaN.
    """
    final_value = traj[-1]
    scaled_distance = (np.linalg.norm(traj - final_value, axis=1)) / np.linalg.norm(traj, axis=1)
    final_index = np.where(scaled_distance < eps)[0]
    if len(final_index) == 0:
        return np.nan
    else:
        return final_index[0]



## Initialize the model
np.random.seed(seed)
eq = GaussianLotkaVolterra(n, random_state=seed, tolerance=1e-10)

## Sample the interaction matrix from a modified normal distribution
eq.A = 1.5 * normal_generalized(size=(n, n), beta=1) - d * np.identity(n)
eq.r = normal_generalized(size=n)

## Random initial conditions
ic = 0.5 * np.random.uniform(size=n)

## Numerically integrate the model
t, y = eq.integrate(tmax, ic)

tind = compute_settling_time(y)
print(f"Settling time: {t[tind]:.2f}")
nonzero_remaining = np.sum(y[-1] > 1e-12) # Precision floor
print(f"{nonzero_remaining} species out of {n} remain")

Settling time: 585.71
146 species out of 300 remain

Hierarchical structure in the Lotka-Volterra model

Trophic interactions can produce extremal values in the Lotka-Volterra model. We consider sampling a family of matrices of the form

$$ A = P^\top (A^{(0)} - d\, I) P + \epsilon \, E, $$

• As before,$A^{(0)}_{ij}, r_i \sim \mathcal{N}(0, 1)$

• $P \in \mathbb{R}^{N \times N}$ is a low-rank assigment matrix describing functional redundancy or trophic levels among species.

• $E \in \mathbb{R}^{N \times N}$ is a matrix of random noise that prevents exact redundancy.

• $\epsilon \ll 1$ is a small constant that controls the degree to which redundant species deviate from exact redundancy.

## Loads a figure from the resources folder
Image("./resources/fig_model.jpg", width=1200)

For an analysis of equilibrium properties of a Lotka-Volterra model with hierarchical structure like this, see Tikhonov (2017).

m_trophs = 3
n_species = 6
tmax = 1e8

interactions_trophic = np.array([
    [0, 0.5, 0],
    [-0.5, 0, 0.3],
    [0, -0.8, 0],
])*50

p_mat = np.array(
    [[1, 0, 0, 0, 0, 0],
     [0, 1, 1, 0, 0, 0],
     [0, 0, 0, 1, 1, 1]]
)


eq = RandomLotkaVolterra(n_species)
eq.r = np.array([0.8, 0.9, 1]) @ p_mat
eq.A = p_mat.T @ (interactions_trophic - 5*np.identity(m_trophs)) @ p_mat #+ 1e-6 * np.random.normal(size=(n_species, n_species))

ic = np.ones(n_species) * 0.1
t, y = eq.integrate(tmax, ic)

plt.semilogx(t, y);
plt.legend(np.arange(n_species), frameon=False);
plt.xlim(1e-3, t[-1]);
plt.xlabel("Time");
plt.ylabel("Species abundance");

Larger ecosystems

• Instead of imposing a specific functional redundancy, we'll sample the assignment matrix $P$ from a family of random low-rank matrices.

• The small parameter $\epsilon$ controls the degree to which redundant species deviate from exact redundancy.

tmax = 1e8
n_val = 200

ic = np.random.uniform(size=n_val)
eq = RandomLotkaVolterra(n_species=200, eps=1e-3)
t, y = eq.integrate(tmax, ic)
plt.semilogx(t, y);
plt.xlim(1e-3, None)

(0.001, np.float64(145205.60974303592))

tmax = 1e8

n_val = 200
ic = np.random.uniform(size=n_val)
settling_times = []
eps_vals = np.logspace(-5, 0, 10)
for eps in eps_vals:
    print(f"eps: {eps}")
    eq = RandomLotkaVolterra(n_val, random_state=0, eps=eps,  sigma=2.0, d=4.5, kfrac=0.2, connectivity=0.05)
    t, y = eq.integrate(tmax, ic)
    tind = compute_settling_time(y)
    settling_times.append(t[tind])


plt.loglog(eps_vals, settling_times, '.')
plt.xlabel("Distance to Exact Functional Redunancy (epsilon)")
plt.ylabel("Settling Time")

eps: 1e-05
eps: 3.5938136638046256e-05
eps: 0.0001291549665014884
eps: 0.0004641588833612782
eps: 0.0016681005372000592
eps: 0.005994842503189409
eps: 0.021544346900318846
eps: 0.07742636826811278
eps: 0.2782559402207126
eps: 1.0

Text(0, 0.5, 'Settling Time')

The Lotka-Volterra model as analogue linear regression

We consider random ecosystems given by the generalized Lotka-Volterra equation,

$$ \frac{dN_i}{dt} = N_i \left( r_i + \sum_{j=1}^N A_{ij} N_j \right) $$

where $N_i$ is the population of species $i$, $r_i$ is the intrinsic growth rate of species $i$, and $A_{ij}$ is the interaction coefficient between species $i$ and $j$. The steady-state solutions of this equation has the form

• Without the outer term, the Lotka-Volterra model is a linear time-invariant system.

• The solution to the Lotka-Volterra model is the solution to the linear system $A \mathbf{N}^* = -\mathbf{r}$.

• The outer term adds a positivity constraint. In practice, this leads to sparsity (exclusion) in the steady-state solution $\mathbf{N}^*$.

eq = RandomLotkaVolterra(n_val, random_state=3, eps=1e0, sigma=0.5, d=4.5, kfrac=0.0, connectivity=0.05)
ic = np.random.uniform(size=n_val)
t, y = eq.integrate(tmax, ic, tol=1e-10)
plt.plot(y);
plt.xlim(1e-3, None)
plt.xlabel("Time")
plt.ylabel("Species abundance")
plt.show()

## Perform unconstrained linear regression
plt.figure(figsize=(6, 3))
plt.plot(np.linalg.solve(eq.A, -eq.r), y[-1], '.')
plt.xlabel("Linear regression solution")
plt.ylabel("Steady-state of Lotka Volterra model")
plt.title("Steady-state solution vs. final time-series value")
plt.show()

## Perform constrained linear regression
from scipy.optimize import nnls, lsq_linear
sol_nnls = lsq_linear(eq.A, -eq.r, bounds=(0, np.inf), tol=1e-10).x
plt.figure(figsize=(4, 4))
plt.plot(sol_nnls, y[-1], '.')
plt.xlabel("Constrained inear regression solution")
plt.ylabel("Steady-state of Lotka Volterra model")
from scipy.stats import spearmanr
corr = spearmanr(sol_nnls, y[-1])
plt.title(f"Correlation: {corr.correlation:.2f}")

Text(0.5, 1.0, 'Correlation: 0.81')

What sets the hardness of linear regression?

We can write a linear regression system of the form $A \mathbf{N} = -\mathbf{r}$ in terms of column vectors of $A$ and $\mathbf{r}$.

If $$ A = \begin{bmatrix} \mathbf{a}_1 & \mathbf{a}_2 & \cdots & \mathbf{a}_N \end{bmatrix} $$

Then we seek a solution to

$$ \mathbf{a}_1 N_1 + \mathbf{a}_2 N_2 + \cdots + \mathbf{a}_N N_N = -\mathbf{r} $$

We write the objective function as

$$ \mathcal{L} = \frac{1}{2} \left( \mathbf{a}_1 N_1 + \mathbf{a}_2 N_2 + \cdots + \mathbf{a}_N N_N + \mathbf{r} \right)^2 $$

If we use gradient descent to solve this system, we can write the update rule as

$$ N_i \leftarrow N_i - \eta \left( \sum_{j=1}^N a_{ij} N_j + r_i \right) $$

Condition: two row-vector products are very

Functional reduncancy is ill conditioning

• Unlike $\epsilon$, which is a property of our particular model for generating interaction matrices $A$, the condition number is a property we can measure from any interaction $A$.

• Two nearly-parallel rows of $A$ is equivalent to two species being nearly-redundant (interchangeable in the community).

tmax = 1e8

n_val = 200
ic = np.random.uniform(size=n_val)
settling_times = []
eps_vals = np.logspace(-5, 0, 10)
condition_numbers = []
for eps in eps_vals:
    print(f"eps: {eps}")
    eq = RandomLotkaVolterra(n_val, random_state=0, eps=eps,  sigma=2.0, d=4.5, kfrac=0.2, connectivity=0.05)
    t, y = eq.integrate(tmax, ic)
    tind = compute_settling_time(y)
    settling_times.append(t[tind])
    condition_numbers.append(np.linalg.cond(eq.A))


plt.figure(figsize=(5, 5))
plt.loglog(eps_vals, settling_times, '.')
plt.xlabel("Distance to Exact Functional Redunancy")
plt.ylabel("Settling Time")

plt.figure(figsize=(5, 5))
plt.loglog(eps_vals, condition_numbers, '.k')
plt.xlabel("Distance to Exact Functional Redunancy")
plt.ylabel("Condition Number")

Text(0, 0.5, 'Condition Number')

tmax = 1e8

n_val = 200
ic = np.random.uniform(size=n_val)
settling_times = []
eps_vals = np.logspace(-5, 0, 10)
condition_numbers = []
for i in range(10):
    progress_bar(i, 20)
    eq = RandomLotkaVolterra(n_species=n_val, random_state=i, eps=np.random.exponential(scale=1e-3), sigma=2.0, d=4.5, kfrac=0.2, connectivity=0.05)
    # eq = GaussianLotkaVolterra(n_val, random_state=seed)
    # eq.A = 1.5 * normal_generalized(size=(n_val, n_val), beta=np.random.exponential(scale=0.5)) - d * np.identity(n_val)
    # eq.r = normal_generalized(size=n_val)
    t, y = eq.integrate(tmax, ic)
    tind = compute_settling_time(y)
    settling_times.append(t[tind])
    condition_numbers.append(np.linalg.cond(eq.A))

plt.figure(figsize=(5, 5))
plt.loglog(condition_numbers, settling_times, '.')
plt.xlabel("Condition Number")
plt.ylabel("Settling Time")

[#########           ] 

Text(0, 0.5, 'Settling Time')

Krylov bound on linear solvers

For an iterative linear solver, the convergence time is given by

$$ \tau \propto \log\bigg( \frac{\kappa(A) - 1}{\kappa(A) + 1} \bigg)^{-1} $$

See Trefethen and Bau (1997).

## Loads a figure from the resources folder
Image("./resources/fig_overview.jpg", width=1200)

How does the settling time depend on the number of species?

n_vals = np.unique(5 * np.logspace(0, 2, 10).astype(int)) # Range of species numbers to consider
tmax = 1e6 # Maximum integration time


## Loop over species numbers and compute the settling time for each
all_t = []
all_y = []
settling_times = []
for i, n_val in enumerate(n_vals):
    progress_bar(i, len(n_vals))

    eq = RandomLotkaVolterra(
        n_species=n_val, 
        random_state=0, 
        eps=1e-3, 
        sigma=2.0, 
        d=6.5, 
        kfrac=0.2, 
        connectivity=0.05
    )
    ic = np.random.uniform(size=eq.n)
    t, y = eq.integrate(tmax, ic)
    all_t.append(t.copy())
    all_y.append(y.copy())
    
    tind = compute_settling_time(y)
    if not np.isnan(tind):
        settling_times.append(t[tind])
    else:
        settling_times.append(np.nan)

[#################   ] 

plt.loglog(n_vals, settling_times, '.k')
plt.xlabel("Number of species")
plt.ylabel("Settling time")
# plt.xlim(None, 1000)
plt.ylim(1e1, None)

(10.0, np.float64(1265833.0390973415))

Supertransients arise in hard optimization problems

• Chimeras in the Kuramoto model
• Pipe turbulence
• Couple map lattices
• Continuous time KSAT solvers

Dimensionality reduction as preconditioning

• What if we project the dynamics into a lower-dimensional space?

• We start from our dynamics in the full $N$-dimensional space

$$ \dot{\mathbf{N}} = \mathbf{f}(\mathbf{N}) $$

• We can instead study the projected by finding a projection operator $P \in \mathbb{R}^{N \times k}$, $k \ll N$, that projects the dynamics into a lower-dimensional space of dimension $k$

  $$ \mathbf{x} \equiv P \mathbf{N}, \quad \dot{\mathbf{x}} = P \dot{\mathbf{N}} $$

  where $\mathbf{x}(t) \in \mathbb{R}^k$ is a lower-dimensional state vector.

## Start by simulating a full-dimensional state vector

n_val = 200
eq = RandomLotkaVolterra(100, random_state=0)
ic = np.random.uniform(size=eq.n)
t, y = eq.integrate(tmax, ic)

plt.figure()
plt.semilogx(t, y);
plt.xlabel("Time")
plt.ylabel("Species abundance")
plt.xlim(1e-2, None)

(0.01, np.float64(794328234.7242821))

How to choose a good projection $P$?

• Since ill-conditioning leads to timescale separation, we seek a projection $P$ satisfying the following properties:

$$ \kappa(P A) \ll \kappa(A) $$

• In numerical analysis, this is known as a preconditioner. If we have a stiff or ill-conditioned linear system, we can precondition it with an appropriate choice of $P$.

• In PDEs, spectral methods can be thought of as a preconditioner. Some problems are ill-conditioned in the physical space, well-conditioned in the Fourier space.

print(f"The condition number of the full system is {np.linalg.cond(eq.A)}")

The condition number of the full system is 19686952.699464306

Singular value decomposition

• Every matrix $A$ can be decomposed into a product of three matrices:

$$ A = U \Sigma V^\top $$

• where $U \in \mathbb{R}^{N \times N}$ and $V \in \mathbb{R}^{N \times N}$ are orthogonal matrices (rotation matrices), and $\Sigma \in \mathbb{R}^{N \times N}$ is a diagonal matrix of positive singular values (in decreasing order).

• We view $U$ as a rotation matrix that rotates our problem into more favorable coordinates.

U, S, V = np.linalg.svd(eq.A)

plt.figure()
plt.plot(S);
plt.xlabel("Rank")
plt.ylabel("Singular value")
plt.title("Singular values of the interaction matrix")
plt.semilogy();

How do we interpret the orthogonal matrices $U$ and $V$?

• The right singular vectors $V_i$ isolate modes with dynamical timescales $\tau_i \sim 1 / \sigma_i$.

## Loads a figure from the resources folder
Image("./resources/fig_svd.jpg", width=1200)

U, S, V = np.linalg.svd(eq.A)
P_fast, P_slow = V[:10, :], V[-10:, :]

yp_fast = P_fast @ y.T # First k modes
yp_slow = P_slow @ y.T # Last k modes


plt.figure()
plt.semilogx(t, y);
plt.xlabel("Time")
plt.ylabel("Species amplitude")
plt.xlim(1e-2, None)


plt.figure()
plt.semilogx(t, yp_fast.T);
plt.xlabel("Time")
plt.ylabel("Fast mode amplitude")
plt.xlim(1e-2, None)

plt.figure()
plt.semilogx(t, yp_slow.T);
plt.xlabel("Time")
plt.ylabel("Slow mode amplitude")
plt.xlim(1e-2, None)
plt.xlabel("Time")
plt.ylabel("Slow mode amplitude")
plt.xlim(1e-2, None)

(0.01, np.float64(2542579.412082511))

Preconditioning: We can compute the condition number of the dynamics projected into the subspaces spanned by the fast and slow modes.

$$ A_\text{slow} = P_\text{slow}^\top A P_\text{slow} $$

For a linear dynamical system, this would also correspond to the Jacobian of the projected system.

print("Original condition number: ", np.linalg.cond(eq.A))

print("Fast modes condition number: ", np.linalg.cond(P_fast @ eq.A @ P_fast.T))
print("Slow modes condition number: ", np.linalg.cond(P_slow @ eq.A @ P_slow.T))

Original condition number:  4429801842.761199
Fast modes condition number:  1.5851693011189856
Slow modes condition number:  153.51209647795702

Why does this work? The condition number alternatively can be expressed as the ratio of the largest and smallest singular values of the matrix.

$$ \kappa(A) = \frac{\sigma_{\text{max}}(A)}{\sigma_{\text{min}}(A)} $$

How general is this? Check out Weyl's inequalities discussed in Thibeault et al. 2024.

When would ill-conditioning emerge in a random ecosystem?

• We revisit the case where the dynamics are not ill-conditioned, the Gaussian Lotka-Volterra model.

eq = GaussianLotkaVolterra(
    n=100, 
    d=8.5,
    random_state=0,
)
tmax = 1000
ic = np.random.uniform(size=eq.n)

t, y = eq.integrate(tmax, ic)

plt.figure()
plt.semilogx(t, y);
plt.xlabel("Time")
plt.ylabel("Species abundance")
plt.xlim(1e-2, t[-1])


print(f"The number of surviving species is {np.sum(y[-1] > 1e-2)} out of {eq.n}")

The number of surviving species is 49 out of 100

• We now perform a simple evolutionary algorithm to evolve the ecosystem to have more surviving species.

• We start with a random interaction matrix, simulate the ecosystem, and find the surviving species.

• If a species is excluded from the ecosystem, we replace it in the next generation with a new species with interactions averaged from two randomly chosen surviving species.

• If $i$ is the index of the extinct species, $j$ and $k$ are the indices of the two surviving species, and $a_{ij}$ is the interaction between species $i$ and $j$, we update the interaction matrix as

  $$ a_{ik} = 0.05 a_{ij} + 0.95 a_{ik} $$

• We repeat this process for 100 generations.

all_fitnesses = []
all_condition_numbers = []
for j in range(100):

    ## Simulate the ecosystem with the current interaction matrix
    t, y = eq.integrate(tmax, ic)
    sol = y[-1]
    
    ## Find the surviving species
    where_survive = (sol > 1 / eq.n)
    surviving_inds = np.where(where_survive)[0] 
    extinct_inds = np.setdiff1d(np.arange(eq.n), surviving_inds)


    ## Recombination: replace extinct species with new species averaged from surviving species
    new_matrix = eq.A.copy()
    for i in extinct_inds:
        parent_ind1, parent_ind2 = np.random.choice(surviving_inds, size=2, replace=False)
        interaction_ind = np.random.choice(np.setdiff1d(np.arange(eq.n), [parent_ind1, parent_ind2, i]))
        new_matrix[i] = eq.A[i]
        new_matrix[i, interaction_ind] = 0.5 * eq.A[parent_ind1, interaction_ind] + 0.5 * eq.A[parent_ind2, interaction_ind]

    ## Update the interaction matrix
    eq.A = new_matrix.copy()
    fitness = np.sum(where_survive)
    print(f"Iteration {j}: The number of surviving species is {fitness} out of {eq.n}")

    ## Save the fitness and condition number for this generation
    all_fitnesses.append(np.sum(where_survive))
    all_condition_numbers.append(np.linalg.cond(eq.A[surviving_inds][:, surviving_inds]))

Iteration 0: The number of surviving species is 49 out of 100
Iteration 1: The number of surviving species is 50 out of 100
Iteration 2: The number of surviving species is 52 out of 100
Iteration 3: The number of surviving species is 53 out of 100
Iteration 4: The number of surviving species is 53 out of 100
Iteration 5: The number of surviving species is 54 out of 100
Iteration 6: The number of surviving species is 55 out of 100
Iteration 7: The number of surviving species is 55 out of 100
Iteration 8: The number of surviving species is 59 out of 100
Iteration 9: The number of surviving species is 59 out of 100
Iteration 10: The number of surviving species is 57 out of 100
Iteration 11: The number of surviving species is 56 out of 100
Iteration 12: The number of surviving species is 58 out of 100
Iteration 13: The number of surviving species is 59 out of 100
Iteration 14: The number of surviving species is 57 out of 100
Iteration 15: The number of surviving species is 59 out of 100
Iteration 16: The number of surviving species is 60 out of 100
Iteration 17: The number of surviving species is 60 out of 100
Iteration 18: The number of surviving species is 60 out of 100
Iteration 19: The number of surviving species is 60 out of 100
Iteration 20: The number of surviving species is 60 out of 100
Iteration 21: The number of surviving species is 60 out of 100
Iteration 22: The number of surviving species is 61 out of 100
Iteration 23: The number of surviving species is 61 out of 100
Iteration 24: The number of surviving species is 61 out of 100
Iteration 25: The number of surviving species is 61 out of 100
Iteration 26: The number of surviving species is 62 out of 100
Iteration 27: The number of surviving species is 62 out of 100
Iteration 28: The number of surviving species is 62 out of 100
Iteration 29: The number of surviving species is 62 out of 100
Iteration 30: The number of surviving species is 62 out of 100
Iteration 31: The number of surviving species is 62 out of 100
Iteration 32: The number of surviving species is 58 out of 100
Iteration 33: The number of surviving species is 58 out of 100
Iteration 34: The number of surviving species is 58 out of 100
Iteration 35: The number of surviving species is 63 out of 100
Iteration 36: The number of surviving species is 64 out of 100
Iteration 37: The number of surviving species is 68 out of 100
Iteration 38: The number of surviving species is 68 out of 100
Iteration 39: The number of surviving species is 68 out of 100
Iteration 40: The number of surviving species is 68 out of 100
Iteration 41: The number of surviving species is 69 out of 100
Iteration 42: The number of surviving species is 70 out of 100
Iteration 43: The number of surviving species is 70 out of 100
Iteration 44: The number of surviving species is 70 out of 100
Iteration 45: The number of surviving species is 70 out of 100
Iteration 46: The number of surviving species is 70 out of 100
Iteration 47: The number of surviving species is 58 out of 100
Iteration 48: The number of surviving species is 64 out of 100
Iteration 49: The number of surviving species is 63 out of 100
Iteration 50: The number of surviving species is 69 out of 100
Iteration 51: The number of surviving species is 69 out of 100
Iteration 52: The number of surviving species is 67 out of 100
Iteration 53: The number of surviving species is 66 out of 100
Iteration 54: The number of surviving species is 67 out of 100
Iteration 55: The number of surviving species is 70 out of 100
Iteration 56: The number of surviving species is 70 out of 100
Iteration 57: The number of surviving species is 70 out of 100
Iteration 58: The number of surviving species is 70 out of 100
Iteration 59: The number of surviving species is 71 out of 100
Iteration 60: The number of surviving species is 71 out of 100
Iteration 61: The number of surviving species is 71 out of 100
Iteration 62: The number of surviving species is 72 out of 100
Iteration 63: The number of surviving species is 72 out of 100
Iteration 64: The number of surviving species is 72 out of 100
Iteration 65: The number of surviving species is 72 out of 100
Iteration 66: The number of surviving species is 73 out of 100
Iteration 67: The number of surviving species is 74 out of 100
Iteration 68: The number of surviving species is 75 out of 100
Iteration 69: The number of surviving species is 75 out of 100
Iteration 70: The number of surviving species is 75 out of 100
Iteration 71: The number of surviving species is 76 out of 100
Iteration 72: The number of surviving species is 77 out of 100
Iteration 73: The number of surviving species is 78 out of 100
Iteration 74: The number of surviving species is 79 out of 100
Iteration 75: The number of surviving species is 80 out of 100
Iteration 76: The number of surviving species is 80 out of 100
Iteration 77: The number of surviving species is 81 out of 100
Iteration 78: The number of surviving species is 81 out of 100
Iteration 79: The number of surviving species is 81 out of 100
Iteration 80: The number of surviving species is 81 out of 100
Iteration 81: The number of surviving species is 80 out of 100
Iteration 82: The number of surviving species is 80 out of 100
Iteration 83: The number of surviving species is 81 out of 100
Iteration 84: The number of surviving species is 81 out of 100
Iteration 85: The number of surviving species is 81 out of 100
Iteration 86: The number of surviving species is 81 out of 100
Iteration 87: The number of surviving species is 81 out of 100
Iteration 88: The number of surviving species is 81 out of 100
Iteration 89: The number of surviving species is 82 out of 100
Iteration 90: The number of surviving species is 82 out of 100
Iteration 91: The number of surviving species is 84 out of 100
Iteration 92: The number of surviving species is 85 out of 100
Iteration 93: The number of surviving species is 85 out of 100
Iteration 94: The number of surviving species is 85 out of 100
Iteration 95: The number of surviving species is 85 out of 100
Iteration 96: The number of surviving species is 85 out of 100
Iteration 97: The number of surviving species is 85 out of 100
Iteration 98: The number of surviving species is 85 out of 100
Iteration 99: The number of surviving species is 85 out of 100

plt.plot(all_fitnesses);

plt.semilogy(all_condition_numbers)
plt.xlabel("Generation")
plt.ylabel("Condition number")

Text(0, 0.5, 'Condition number')

## Loads a figure from the resources folder
Image("./resources/fig_evolve.jpg", width=1200)

How do transients depend on the initial conditions?

tmax = 1e8 ## Maximum integration time

eq = RandomLotkaVolterra(
    n_species=n_val, # Number of species
    sigma=2.0, # Amplitude of interactions
    d=4.5, # Density-limitation strength
    kfrac=0.5, # Fraction redundant interactions
    eps=0.00005, # Functional redundancy
    random_state=0, # Random seed
    connectivity=0.05 # Connectivity of network
)

ic = 2*np.random.uniform(size=eq.n)
tlim = (0, 50000)
t, y = eq.integrate(tmax, ic)
plt.figure()
plt.semilogx(t, y, color="r", linewidth=0.8);

ic += 1e-3
t, y = eq.integrate(tmax, ic)
plt.semilogx(t, y, color="k", linewidth=0.8);

plt.xlabel("Time")
plt.ylabel("Species abundance")
plt.xlim(1e-2, t[-1])

# settling_time = compute_settling_time(y)
# print(f"Settling time: {t[settling_time]:.2f}")

n1_val = np.linspace(0, 1, 100)
n1_val = np.linspace(0.2, 0.3, 100)
n1_val = np.linspace(0.2, 0.21, 100)
icp = ic.copy()
all_sol = []
for i, n1 in enumerate(n1_val):
    progress_bar(i, len(n1_val))
    icp[0] = n1
    t, y = eq.integrate(tmax, icp)
    all_sol.append([y, t])



    

U, S, V = np.linalg.svd(eq.A)
# P = U[:, :10].T
# P1 = U[:, -10:].T
P_fast, P_slow = V[:10, :], V[-10:, :]
# P_slow = P_slow.T
def instantaneous_fli(y, t, dt, jac):
    
    # P2 = P_slow @ jac(t, y) @ P_slow.T

    # # ## Finite difference approximation of Jacobian
    eps = 1e-8
    jac_fd = np.zeros((P_slow.shape[0], P_slow.shape[0]))
    for i in range(jac_fd.shape[0]):
        pt0, pt1 = y.copy(), y.copy()
        pt1[i] += eps
        jac_fd[:, i] = (P_slow @ (eq(0, pt1) - eq(0, pt0))) / (P_slow @ (pt1 - pt0))
        # jac_fd[:, i] = (P_fast @ (eq(0, pt1) - eq(0, pt0))) / (P_fast @ (pt1 - pt0))
    jac_fd = np.array(jac_fd)
    P2 = jac_fd

    # Finite‐difference reduction:
    
    # eps = 1e-8
    # n_slow = P_slow.shape[1]   # should be 10
    # J_reduced_fd = np.zeros((n_slow, n_slow))
    # f0 = eq(0, y)  # baseline 100‐vector
    # for j in range(n_slow):
    #     p_j = P_slow[:, j]              # 100‐vector (slow direction j)
    #     y_perturbed = y + eps * p_j     # perturb along the jth slow basis
    #     f1 = eq(0, y_perturbed)         # 100‐vector
    #     # directional derivative of f in direction p_j:
    #     df_dp = (f1 - f0) / eps         # 100‐vector ≈ J_full @ p_j
    #     # project onto slow‐coordinates on the left:
    #     J_reduced_fd[:, j] = P_slow.T @ df_dp   # 10‐vector
    # P2 = J_reduced_fd

    # Let U[k, :] be the k-th slow‐mode represented as a 100‐vector.
    # U = P_slow.T       # shape (10, 100)
    # U = P_slow.T
    # J_full = jac(t, y)
    # # print(J_full.shape, U.shape)
    # # Preallocate the “scaled‐projection” matrix:
    # #   jac_fd[k, j] = (U[k, :] · (J_full[:, j])) / U[k, j]
    # n_modes, dim = U.shape
    # jac_fd  = np.zeros((n_modes, dim))

    # numerator = U @ J_full  # (10, 100)
    # # We only need the first n_modes columns (j=0..n_modes-1) to match the original FD over 10 coordinates.
    # # Denominator = U[:, :n_modes] (shape: 10 × 10)
    # num_slice = numerator[:, :n_modes]  # (10, 10)
    # den_slice = U[:, :n_modes]       # (10, 10)
    # # Compute jac_fd[k, j] = (u_k^T · (J_full · e_j)) / u_{k,j} for j=0..n_modes-1
    # jac_fd = num_slice / den_slice  # shape: (10, 10)
    # P2 = jac_fd

    

    Q = np.linalg.inv(np.eye(P2.shape[0]) - dt * P2)
    eigenvalues = np.linalg.eigvalsh(Q.T @ Q)
    fli_val = np.sqrt(np.max(eigenvalues))
    return fli_val


def find_fli(sol, tvals, jac):
    dtvals = np.diff(tvals)
    dtvals = np.hstack([dtvals, dtvals[-1]])
    fli = -np.inf
    for yval, dt, tval in zip(sol, dtvals, tvals):
        fli = max(fli, instantaneous_fli(yval, tval, dt, jac))
    return fli

all_fli = []
for sol, tvals in all_sol:
    all_fli.append(find_fli(sol, tvals, eq.jac))

plt.plot(n1_val, all_fli);

# n_res = 120
# n_res = 100
# n_res = 10
n_res = 10
# n_res = 10
# n_res = 25
# n_res = 10
# n_res = 5
# x = np.linspace(0.2, 0.202, n_res)
# y = np.linspace(0.2, 0.202, n_res)
x = np.linspace(0.2, 0.3, n_res)
y = np.linspace(0.2, 0.3, n_res)
x = np.linspace(0.0, 1.0, n_res)
y = np.linspace(0.0, 1.0, n_res)
xx, yy = np.meshgrid(x, y)
xx, yy = xx.ravel(), yy.ravel()
pts = np.vstack([xx, yy]).T

## Fixed disorder in initial conditions
np.random.seed(0)
# ic0 = np.random.uniform(size=eq.n)
icax = np.random.uniform(size=eq.n) * np.max(sol_final)
icbx = np.random.uniform(size=eq.n) * np.max(sol_final)
icay = np.random.uniform(size=eq.n) * np.max(sol_final)
icby = np.random.uniform(size=eq.n) * np.max(sol_final)



all_fli = []
for i, ic_val in enumerate(pts):
    progress_bar(i, len(pts))
    ic = bilinear_interpolation(icax, icbx, icay, icby, ic_val[0], ic_val[1])
    ic[ic < 0] = 0.0
    # print(ic[0])
    t, y = eq.integrate(tmax, ic, tolerance=1e-10)
    fli = find_fli(y, t, eq.jac)
    all_fli.append(fli)
    

ms = 1200 * (20 / n_res)**2
plt.figure(figsize=(8,8))
all_fli = np.array(all_fli).astype(float)
vmin = np.percentile(all_fli, 5)
vmax = np.percentile(all_fli, 95)
plt.scatter(pts[:, 0], pts[:, 1], c=all_fli, s=ms, vmin=vmin, vmax=vmax)
# plt.colorbar()
plt.axis("off")

## Loads a figure from the resources folder
Image("./resources/fig_chaos.jpg", width=1200)

The ill-conditioned Lotka-Volterra model

We consider random ecosystems given by the generalized Lotka-Volterra equation,

$$ \frac{dN_i}{dt} = N_i \left( r_i + \sum_{j=1}^N A_{ij} N_j \right) $$

where $N_i$ is the population of species $i$, $r_i$ is the intrinsic growth rate of species $i$, and $A_{ij}$ is the interaction coefficient between species $i$ and $j$. The steady-state solutions of this equation has the form

$$ -A \mathbf{N}^* = \mathbf{r} $$

where $N_i \geq 0$ for all $i$.

In this notebook, we explore the behavior of this model when the interaction matrix $A_{ij}$, and the growth rates $r_i$ are drawn from random distributions. Specifically, we consider the case where $r_i \sim \mathcal{N}(0,1)$, and the matrix $A$ has the form

$$ A = P^\top (Q - d\, I) P + \epsilon E $$

where $Q_{ij} \sim \mathcal{N}(0,1)$, $E_{ij} \sim \mathcal{N}(0,1)$, $P$ is a low-rank matrix imposing functional redundancy, $d$ is a constant density-limitation, and $\epsilon \ll 1$ is a small constant.

from scipy.integrate import solve_ivp

## Load local functions
import sys
import base
import utils
from base import RandomLotkaVolterra
from utils import levenshtein


tlim = (0, 20000) # Integration interval
n_val = 200 # Number of species

## Initialize the model
eq = RandomLotkaVolterra(n_val, sigma=2.0, d=4.5, kfrac=0.2, eps=0.001, random_state=0, connectivity=0.01)

## Initial conditions
ic = np.random.uniform(size=eq.n)

## Numerical integration
fsol = solve_ivp(eq, tlim, ic[:n_val], jac=eq.jac, **eq.integrator_args)


## Check that the found solution is stable
sol_final = fsol.y[:, -1]
jac_final = eq.jac(0, sol_final)
eigs = np.linalg.eigvals(jac_final)
print(f"Numerical stability observed: {np.all(np.real(eigs) < 0)}")

## Plot the solution
plt.figure(figsize=(7, 4))
plt.semilogx(fsol.t, fsol.y.T, color="k", lw=1.5, alpha=0.3);
plt.xlabel("Time")
plt.ylabel("Species abundance")
plt.xlim(1e-2, np.max(fsol.t))
plt.ylim(0, None)

eigvals = np.linalg.eigvals(jac_final)
plt.figure(figsize=(4, 4))
plt.plot(np.real(eigvals), np.imag(eigvals), '.k')
plt.xlabel("Real part")
plt.ylabel("Imaginary part")
plt.axis("equal")
plt.title("Eigenvalues of Jacobian")

Numerical stability observed: True

Text(0.5, 1.0, 'Eigenvalues of Jacobian')